Journal  of  Power  Sources  194  (2009)  190-198 


ELSEVIER 


Contents  lists  available  at  ScienceDirect 

Journal  of  Power  Sources 

journal  homepage:  www.elsevier.com/locate/jpowsour 


pri 

SllbiidtS 


A  droplet  size  dependent  multiphase  mixture  model  for  two 
phase  flow  in  PEMFCs 

Guangli  Hea,  Yohtaro  Yamazaki3*,  Abuliti  Abudulab 

a  Department  of  Innovative  and  Engineered  Materials,  Interdisciplinary  Graduate  School  of  Science  and  Engineering,  Tokyo  Institute  of  Technology,  Tokyo,  Japan 
b  New  Energy  Technology  Research  Division,  Aomori  Industrial  Research  Center,  Aomori,  Japan 


ARTICLE  INFO 


ABSTRACT 


Article  history: 

Received  16  January  2009 

Received  in  revised  form  24  March  2009 

Accepted  11  May  2009 

Available  online  20  May  2009 


Keywords: 

Droplet 

Multiphase  mixture  model 
Two  phase  flow 
PEMFC 
Fuel  cell 


A  droplet  size  dependent  multiphase  mixture  model  is  developed  in  this  paper,  and  the  droplet  size  in 
the  gas  channel  can  be  considered  as  a  parameter  in  this  multiphase  mixture  model,  which  includes  the 
effect  of  gas  diffusion  layer  (GDL)  properties  and  the  gas  drag  function  and  cannot  be  considered  in  the 
commonly  used  multiphase  mixture  model  in  the  references.  The  three-dimensional  two  phase  and  non- 
isothermal  simulation  of  the  PEMFCs  with  a  straight  flow  field  is  performed.  The  effect  of  droplet  size  on 
the  liquid  remove,  the  effect  of  liquid  water  on  the  heat  transfer  and  the  effect  of  gas  flow  pattern  on  the 
heat  and  mass  transfer  are  mainly  investigated.  The  simulation  results  show  that  the  large  droplet  is  hard 
to  be  dragged  by  the  gas,  so  it  produces  large  water  saturation.  The  results  of  the  heat  transfer  show  that 
the  liquid  water  hinders  the  heat  transfer  in  the  GDL  and  catalyst  layer,  so  it  produces  the  large  relative 
high  temperature  area,  and  there  are  large  temperature  difference  and  water  saturation  in  the  PEMFCs 
operated  with  coflow  pattern  compared  with  counter  flow  pattern. 

©  2009  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Water  management  is  very  important  for  the  operation  of 
PEMFCs.  Keeping  the  membrane  with  enough  water  content  and 
removing  the  excessive  liquid  water  effectively  are  the  main  focus 
of  water  management,  which  directly  affect  performance  and  heat 
management  of  the  PEMFCs.  Liquid  water  transport  in  PEMFCs 
occurs  as  follows:  (1)  water  is  produced  in  cathode  catalyst  layer, 
and  liquid  water  transports  within  the  gas  diffusion  layer  (GDL) 
by  capillary-driven  flow.  (2)  Liquid  water  droplets  appear  on  the 
GDL/gas  channel  interface  and  are  removed  by  the  gas  shearing 
function  [1-4].  (3)  Liquid  water  travels  in  the  gas  channel.  Wang 
and  his  coworkers  [3]  have  observed  the  emergence,  growth  and 
detachment  of  liquid  water  droplet  on  the  GDL/gas  channel  inter¬ 
face.  And  they  [5]  have  measured  the  size  of  the  droplet  on  the 
GDL/gas  channel  interface,  and  find  that  liquid  water  can  transport 
through  the  gas  channel  without  interaction  with  channel  wall  at 
high  gas  velocity.  Also,  Kimball  et  al.  [6,7  ]  have  proved  that  the  liquid 
water  flow  through  the  largest  pore  of  the  hydrophobic  gas  diffu¬ 
sion  layer,  in  their  experiments,  it  was  found  that  the  liquid  water 
appeared  at  the  same  position  for  the  different  operating  condition, 
which  shows  the  new  way  for  the  analysis  and  insight  knowing  of 
the  liquid  transportation  in  the  gas  diffusion  layer. 
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By  now,  many  models  have  been  developed  for  simulating  the 
liquid  water  transport  in  PEMFCs.  Those  models  are  based  on  the 
different  theory  and  assumption  in  flow  dynamics.  One  of  the  most 
commonly  used  model  is  the  multiphase  mixture  (M2)  model, 
which  consider  the  liquid  and  the  gas  as  mixture  [8-13]  simply 
and  calculation  cost  effectively.  However,  two  phase  flows  in  PEM¬ 
FCs  are  complex  dynamic  processes.  Gurau  et  al.  [13]  show  that 
the  M2  model  has  a  narrow  range  of  applicability,  which  is  limited 
to  steady-state  flows  without  change  of  phase  and  without  phase 
production  due  to  other  physical  processes.  For  more  complex  situ¬ 
ations,  including  those  commonly  encountered  in  PEMFCs,  the  M2 
model  ceases  to  reflect  the  principles  and  could  lead  to  predictions 
of  unrealistic  velocity  and  scalar  fields.  Furthermore,  the  M2  model 
could  represent  a  tool  unable  to  capture  complex  fuel  cell  phe¬ 
nomena  such  as  water  transfer  and  the  liquid  droplet  in  the  gas 
channel  and  its  effect  on  the  mass  transfer  cannot  be  considered  in 
the  multiphase  mixture  model.  In  the  recent  works,  the  two-fluid 
model  has  obtained  more  attention  for  its  convenience  of  consid¬ 
ering  the  liquid  phase  and  gas  phase  separately,  so  it  can  describe 
more  phenomena  in  the  two  phase  flow  [14-16],  but  there  are  also 
some  disadvantages  for  the  two-fluid  model,  such  as  it  needs  more 
calculation  cost,  and  its  hard  to  converge  for  it  includes  the  liquid 
momentum  equation,  liquid  continuity  equation,  gas  momentum 
equation,  and  gas  continuity  equation. 

In  this  paper,  we  developed  a  droplet  size  dependent  multi¬ 
phase  mixture  model  for  considering  the  two  phase  behavior  in 
PEMFCs.  In  this  model,  the  interacting  effect  between  two  phases  is 
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Nomenclature 

Parameters  and  variables 
a  acceleration 

CD  drag  coefficient 

Cj  mass  transfer  coefficient 

D  diffusion  coefficient  (cm2  s-1 ) 

dp  characteristic  droplet  size  (cm) 

F  Faraday’s  constant 

/drag  drag  function 

M  molar  mass  (kg  mol-1 ) 

Mi  interface  force  of  liquid  (N) 

Mg  interface  force  of  gas  (N) 

n  number  of  electrons 

P  pressure  (Pa) 

q  switch  function 

R  universal  gas  constant  (J  mol-1 1<-1 ) 

Peat  reaction  rate  in  cathode 

Pan  reaction  rate  in  anode 

Re  Reynold  number 

s  water  saturation 

T  temperature  (I<) 

V  velocity  (cm  s-1 ) 

Voc  open  circuit  voltage  (V) 

y  molar  fraction 

Greek  letters 
0  potential  (V) 

s  volume  fraction 

a  surface  tension  (N  cm-1) 

6C  equilibrium  contact  angle  on  diffuser 

A  polymer  water  content  H20/S03 

ad  drag  coefficient  of  water  in  membrane 
r\  overpotential  (V) 

fi  viscosity  (Pas) 

p  density  (kg cm-3) 

r p  relaxation  time 

Subscripts 
an  anode 

C  about  capillary 

cat  cathode 

dr  drift 

sat  saturated 

sol  about  electron 

g  gas  phase 

H2  hydrogen 

i  note  for  species 

k  note  for  phases 

l  liquid  phase 

m  mixture  properties  of  multiphase  mixture 

mem  polymer  phase 

02  oxygen 

p  phase 

q  pahse 

P  relative 

w  water 

wv  water  vapor 


considered  in  a  multiphase  mixture  form,  which  includes  the  effect 
of  the  droplet  size,  the  drag  coefficient,  velocity,  Reynolds  number, 
the  droplet  relaxation  time.  Compared  with  the  two-fluid  model, 
the  present  model  is  calculation  cost  effective  and  it  is  easy  to  be 
performed,  also  it  can  includes  the  effect  of  the  droplet  size  on 


Fig.  1.  Geometry  model  for  the  simulation. 

the  liquid  removal  compared  with  the  mixture  model,  which  is  a 
key  factor  in  two  phase  flow  in  PEMFCs.  However,  the  gas  channel 
wall  contact  angle  can  affect  the  liquid  form  in  the  gas  channel, 
and  the  hydrophobic  gas  channel  wall  likely  produces  droplet,  but 
hydrophilic  channel  wall  likely  produces  water  film.  The  interaction 
between  gas  and  the  liquid  film  is  more  complex  than  that  of  water 
droplet,  that  is  will  be  considered  in  the  future  work.  So,  the  present 
model  is  preferable  for  PEMFCS  with  hydrophobic  gas  channel. 

2.  Model  development 

2.1  Model  assumption 

The  calculated  regions  consist  of  conventional  straight  channels, 
gas  diffusion  electrodes,  catalyst  layers  and  a  membrane,  as  well 
as  current  collectors,  which  are  shown  in  Fig.  1.  Because  PEMFC 
is  operated  in  temperature  of  80  °C,  so  it  is  assumed  that  water  is 
generated  in  cathode  catalyst  layer  as  liquid,  and  the  mss  transfer 
between  liquid  phase  and  gas  phase  are  considered.  Because  it  is  a 
steady  state  model  in  this  work,  so  the  growing  up  process  for  the 
droplet  cannot  be  considered,  and  also,  the  effect  of  gas  channel 
wall  is  not  involved  in  this  model.  In  the  gas  diffusion  layer  and 
catalyst  layer,  the  liquid  driven  force  is  capillary  force,  which  is  the 
function  of  water  saturation,  so  the  other  forces  are  not  considered. 

2.2.  Governing  equations 

The  droplet  size  dependent  multiphase  mixture  model  is  devel¬ 
oped  as  follows: 

The  continuity  equation  for  the  multiphase  mixture  is 

9 

T^ipm)  +  V  •  (pmUm)  =  0  (1) 

Vm  =  TLWPWk  (2) 

Pm 

n 

Pm  =  Jfoikpk  (3) 

k= 1 

The  momentum  equation  for  the  mixture  can  be  expressed  as: 

d  _ 

T^iPmVm)  +  V  •  (PmVmVm)  =  —  Vp  +  V  •  [/Xm( V Vm  +  V Vm )]  +  Pmg 

+  F  +  V  '  ^^a/cP)cC’dr,kC’dr,k  j  (4) 
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Mm  = 


k= l 

udr tk  is  the  drift  velocity  for  secondary  phase  k: 

Air ,k  —  Vk  ~ 

The  relative  velocity  is 

Vpq  =  Vp  —  Vq 

The  mass  fraction  for  any  phase  (/<)  is  defined  as: 
akPk 


Ck  = 


Pm 


(5) 


(6) 


(7) 


(8) 


The  drift  velocity  and  the  relative  velocity  (uqp)  are  connected  by 
the  following  expression: 


Alr,p  —  Vpq  ^  ]ck  Vqk 
k=  1 

rp  (Pp  ~  Pm) 


upq  '■ 


/drag  Pp 


a 


(9) 

(10) 


where  rp  is  the  particle  relaxation  time  [17]. 

Particle  relaxation  time  in  two  phase  flow  means  the  time  in 
which  the  droplet’s  velocity  change  from  zero  to  the  equilibrium 
velocity  with  the  gas  phase.  And 


tn  = 


Ppdp 

18  flq 


(11) 


The  default  drag  function  /drag  is  taken  from  Schiller  and  Neu¬ 
mann  [18]: 


Re  <  0.01 
0.01  <Re<  20 
20  <  Re  <  260 


r  _9  24 

D  2  +  Re 

CD  =  ^[1  +0.1315Re(0-82“005w)] 
Re 

CD=  ^[1  +0.1935Re°-6305] 


(12) 


260  <  Re  <  1 .5  x  103  log10  CD  =  1 .6435 

—  1 .1242w  +  0.1558W2 
w  =  log10  Re 

and  the  acceleration  a  is  of  the  form  [19]: 

->  ->  .  _>  .  _>  dUffi 

a  =  g-  (um  •  V)vm  -  -gp 

From  the  continuity  equation  for  phase  p,  the  volume  fraction  equa¬ 
tion  for  phase  p  can  be  obtained: 

d  n 

-T^(oipPp)  +  V  •  ( OipppVm )  =  —V  •  ((XpPpVdr,p)  +  ^  ](fflqp  —  ^pq) 


(13) 


q=i 


(14) 


In  the  gas  difusion  layer  and  catalyst  layer,  the  liquid  is  driven  by 
capillary  force,  then: 


d{spis) 

dt 

And 

Pc  = 


+  V- 


Ks 3  dpc 

A - 4-  Vs 

/X/  ds 


—  Tw 


g  C0S„^  (1.417(1  -s)~  2.12(1  -  s)2  +  1.263(1  —  s)3 )  0C  <90° 

C If/®)0' 5 


<T  COS#c 


(W 


oy  (1.417s  -  2.12s2  +  1.263s3) 


0C  >  90° 


(15) 


(16) 


The  mass  transfer  rate  rw  between  two  phases  is 

cr  max  Q(1  -  s)Pwv~TPsatMw,H2o]  ,  [-8p,])  (17) 

The  species  mass  conservation  in  phase  p  is 

V  •  (ctgpg~i? gy{)  =  V  •  (D/Vyj)  +  S,-  (18) 

And,  Vg  can  be  obtained  according  to  Eqs.  (6),  (7)  and  (10),  source 
term  for  the  gas  phase  is 

_  Mw,  h2 
*^h2  —  ^  Rz 


So2=- 


M, 


2F 

Fw,02 

4F 


f^cat 


c  MW,H20  d 

^h2o  =  — 2f — Kcat 


(19) 

(20) 
(21) 


In  catalyst  layer  the  relationship  between  species  mass  fraction  on 
the  catalyst  surface  and  the  reaction  sites  is 


(22) 


(23) 


PO I  ly  _y  N  _  MWJ 

^  UP,  surf  yi,cent)r  —  ^  ^an,cat 

The  effective  diffusion  coefficient  is  [20]: 

The  membrane  phase  and  solid  phase  potential  conservation  equa¬ 
tions  and  electrochemical  reaction  rate  in  the  cathode  side  and 
anode  side: 

^  '  (^sol^^sol)  +  Ksol  —  9  (24) 

V  •  (CTmem V 0mem  )  +  ^mem  =  0  (25) 

Ran  =  Jan  (  jyy-  )  ^  (e“anF,'an/RT  -  e-““t^an/RT )  (26) 

Kcat  =ja[  nFricx,RT  ~  e~aacFr)c3^RT)  (27) 

where 

Pan  =  0sol  —  0mem  (28) 

Peat  =  0sol  —  0mem  —  Aoc  (29) 

For  the  consideration  of  the  effect  of  liquid,  the  reaction  rate  calcu¬ 
lated  in  Eqs.  (26)  and  (27)  are  multiplied  by  (1  -s;). 

Membrane  phase  electric  conductivity  [21  ]: 


<7mem  =  yfe(0.514A  -  0.326)V268((1/303)-l/T) 


a«  =  2-522 


X  =  0.043  +  17.18a-  39.85a2  +36a3(a  <  1) 

A.  =  14  +  1 .4(a  -  1  )(a  >  1) 
a=^2s 

Ksat 

Pwv  =  ^H2oK 

logto  Psat  =  -2.1794  +  0.02953(T-  273.17)  -9.1837 


(30) 

(31) 

(32) 

(33) 

(34) 


x  10"5(T  -  273.17)z  +  1.4454  x  10“'(7  -  273.17) 


,-7,- 


Water  diffusion  flux  through  the  membrane  is 
Jwff  =  -£rMh20  D,VX 


(35) 

(36) 
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And  [21  ] 

D,  =/(A)e2416«1/  3Q3)-i/n  (37) 

The  heat  conservation  equation  is 
V  ■  (PmCpm^mT)  =  V  •  (fcejy-VT)  +  S/| 

i=n 

T  =  sr^^iPJj 

pm 
i= 1 

^  (38) 

fee jf  =  £}  OCjki  +  (1  -  £)fes 

i=l 

i=n 

cpm  =  ^  '(XjCpj 

i=l 

5h  is  calculated  as  follows: 

Sjj  =  I  R0hm  +  ^reaction  +  ^  an, cat  +  frphase  (39) 

2.3.  Boundary  conditions  and  parameters 

The  inlet  volume  flow  rate  is  150  cm3  min-1 ,  which  is  converted 
to  mass  flow  rate  by  the  UDF  (user  defined  function  in  Fluent® 
software),  the  outlet  boundary  condition  is  the  pressure  outlet  con¬ 
dition,  the  outlet  pressure  is  equal  to  atmosphere  pressure.  The  inlet 
temperature  is  353  K,  and  the  temperature  of  anode  current  collec¬ 
tor  and  cathode  current  collector  end  walls  boundary  are  the  353  K, 
which  means  the  current  collector  is  ideally  cooled.  The  other  lat¬ 
eral  walls  and  the  end  walls  are  impermeable  for  all  the  species.  The 
operated  potential  is  set  on  the  boundary  of  cathode  current  col¬ 
lector.  While  the  potential  on  the  anode  current  collector  boundary 
is  set  to  be  0.  The  water  saturation  in  cathode  and  anode  inlet  is 
zero.  In  the  two-fluid  model  and  the  droplet  size  dependent  multi¬ 
phase  mixture  model  in  the  present  study,  the  droplet  size  in  the  gas 
channel  is  the  key  variable.  Which  integrate  the  effect  of  the  prop¬ 
erties  of  GDL  and  gas  flow  properties  on  the  formation  of  liquid 
droplet,  including  the  commonly  used  contact  angle,  and  the  pore 
size  characters  of  the  GDL,  as  well  as  the  structure  of  GDL.  There  are 
some  works  which  have  concentrated  on  the  determination  of  the 
droplet  detachment  diameter  by  the  computational  fluid  dynam¬ 
ics  (CFD)  method  [22-24]  and  analytical  method  [14],  but  the  GDL 
structure  cannot  be  considered  in  all  those  models,  which  is  essen¬ 
tial  for  water  flow  through  the  GDL  to  the  gas  channel  [6,7].  So,  the 
experiential  formula  for  estimating  the  droplet  size  is  applied  as 
follows,  which  was  completed  by  Wang  and  co-worker  [5]  for  the 
certain  experimental  condition: 

log(dp)  =  — 2.591og(i ;)  +  Jf  -  1.591og(l  +  5.2i?e“0-63)  (40) 

The  values  of  the  other  parameters  used  in  the  model  are  listed  in 
Table  1. 

2.4.  Mesh  grid  and  solution  technique 

The  geometry  model  is  shown  in  Fig.  1  and  it  is  discretized  into 
550,500  hexahedral  mesh  volumes,  and  to  assure  the  quality  of  the 
grid,  the  size  of  the  grid  in  gas  channel,  gas  diffusion  layer,  catalyst 
layer  and  membrane  are  different.  The  Simplec  algorithm  and  Quick 
difference  scheme  are  applied  for  solving  the  pressure-velocity 
coupled  equations,  and  species  equations.  And  suitable  relax  factors 
are  used  for  momentum,  slip  velocity,  water  saturation,  potential 
and  species.  The  simulation  is  performed  in  Flunet®  software  of 
Ansys  company  with  some  codes  of  UDF  (User  Defined  Function  in 
Fluent)  added  by  us. 


Table  1 

Values  of  the  parameters. 


Physical  properties 

Value 

Faraday’s  constant,  F 

96487  C  moU1 

Permeability  of  gas  diffusion  layer,  I<p 

8  x  10“8  cm2 

Liquid  water  viscosity,  /x; 

3.565  x  10“4 Pas 

Anodic  transfer  coefficient,  aa 

0.5 

Cathodic  transfer  coefficient,  ac 

0.55 

Water  contact  angle  in  diffuser,  0 

120° 

Gas  channel  width 

0.1  cm 

Gas  channel  length 

5  cm 

Gas  channel  height 

0.1  cm 

Thickness  of  current  collector 

0.15  cm 

Anode  GDL  thickness 

0.019  cm 

Cathode  GDL  thickness 

0.019  cm 

Gas  diffusion  layer  void  fraction 

0.7 

Catalyst  layer  thickness 

0.002  cm 

Catalyst  layer  void  fraction 

0.5 

Membrane  thickness  (Nafion®112) 

0.0005  cm 

Cell  inlet  temperature 

353  K 

Outlet  pressure 

0.1  MPa 

Air  and  fuel  inlet  humidified  temperature 

348  K 

Open  circuit  voltage 

0.95  V 

Mass  transfer  rate  between  phases,  cT 

200  s-1 

Gas  constant,  R 

8314J  mol-1  K"1 

Reference  hydrogen  concentration,  HT2ef  1  kmol  m-3 
Reference  oxygen  concentration,  O^ 

1  kmol  m3 

Operation  current  density 

1.0  A  cm-2 

Anode  exchange  current  density,  jjjf 

1.5e8  Am-3 

Cathode  exchange  current  density, 

7000  Am-3 

Thermal  conductivity  of  GDL,  membrane  and  CL,  k 

8  wirr1  k-1 

Thermal  current  collector 

8  WITT1  k-1 

Specific  heat  capacity  of  liquid  water 

4182  J  kg-1  k-1 

Electrical  conductivity  of  GDL  and  CL 

5000 1  ohm-1  itt1 

3.  Results  and  discussion 

3.1.  Validity  of  the  present  work 

In  many  references,  the  model  of  PEMFCs  were  validated  by  the 
comparison  between  the  calculated  performance  and  the  experi¬ 
ment  data,  although  it  is  impossible  to  account  for  all  the  factors 
for  the  performance  in  the  model,  there  are  still  no  more  suitable 
method  for  validating  the  model  results.  In  some  references,  the 
liquid  water  distribution  was  observed  by  optical  method  or  neu¬ 
tron  tomogram,  but  it  is  just  the  qualitative  method.  So,  to  explore 
the  simple  and  precise  model  validated  method  is  also  a  subject 
which  should  be  paid  more  attention  in  the  future  work.  To  show 
the  reasonability  of  the  present  model,  the  experimental  data  and 


Fig.  2.  Comparison  between  the  experimental  data  and  simulation  result. 
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0.0020 


Fig.  3.  Oxygen  mass  fraction  in  the  cathode  side  for  the  base  case. 


the  calculated  result  from  PEMFCs  Module  (one  of  the  add-in  mod¬ 
ule)  in  Fluent®  software  are  both  used.  Fig.  2  shows  the  comparison 
among  the  three  results.  It  can  be  seen  that  the  present  model  is 
more  suitable  than  the  PEMFCs  Module  for  the  high  current  den¬ 
sity.  Because  the  liquid  is  assumed  to  have  the  same  velocity  and 
only  gas  momentum  equation  is  solved  in  the  PEMFCs  Module  in 
Fluent  software. 

3.2.  Species  mass  fraction,  water  saturation,  and  current  density 
in  the  PEMFCs 

Figs.  3-5  show  the  oxygen  mass  fraction,  water  saturation 
and  current  density  distribution  in  the  PEMFCs,  respectively.  As 
expected,  the  oxygen  mass  fraction  decreases  from  inlet  to  outlet, 
from  gas  channel  to  catalyst  layer.  And  the  low  oxygen  mass  fraction 
area  appears  in  the  catalyst  layer  and  GDL  neighboring  to  the  rib  of 
the  current  collector,  where  the  longest  distance  exists  for  the  gas 
convection  and  diffusion. 

The  same  distribution  characters  are  obtained  for  the  water  sat¬ 
uration  shown  in  Fig.  4.  However,  liquid  water  behavior  in  the  gas 
flow  channels  and  GDL  are  very  complicated,  especially  for  the 
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Fig.  5.  Current  density  distribution  in  the  PEMFCs  for  the  base  case,  A  m-2. 


droplets  in  the  gas  flow  channels,  when  the  size  of  the  droplet 
reaches  a  certain  value,  it  may  interact  with  the  gas  channel  walls 
[25,5],  then  the  liquid  droplet  maybe  deform  to  liquid  film  if  the 
channel  wall  is  hydrophilic.  Which  are  difficult  to  be  considered  in 
the  mixture  model  or  two-fluid  model.  And  the  results  obtained  in 
the  two-fluid  model  or  the  mixture  model  for  the  water  saturation 
are  the  volume-averaged  values.  So,  the  water  saturation  obtained 
in  the  GDL  is  more  meaningful  than  that  in  the  gas  channel.  And 
it  can  be  observed  in  Fig.  4  that  the  current  collector  ribs  hinder 
the  liquid  water  removal  so  the  water  saturation  in  the  GDL  and 
catalyst  layer  neighboring  to  the  ribs  is  higher  compared  with  that 
neighboring  the  gas  channel. 

Fig.  5  shows  the  current  density  distribution  in  the  whole 
PEMFCs,  it  can  be  seen  that  the  highest  current  density  appears  in 
the  corner  between  the  gas  channel  and  current  collector  ribs.  Due 
to  the  current  cannot  be  conducted  through  the  channel  so  the 
current  density  is  low  in  the  area  neighboring  the  gas  channel  top 
and  bottom. 


3.3.  Effect  of  water  saturation  on  the  temperature  distribution 


Fig.  6  shows  the  temperature  distribution  in  the  PEMFCs  for 
the  base  case  (the  operation  current  density  is  about  1.0  A  cm-2). 


355.687 

355.545 

355.404 

355.263 

355.121 

354.980 

354.838 

354.697 

354.556 

354.414 

354.273 

354.131 

353.990 

353.848 

353.707 

353.566 

353.424 

353.283 

353.141 

353.000 


Fig.  4.  Water  saturation  in  the  cathode  side  for  the  base  case. 


Fig.  6.  The  temperature  distribution  in  the  PEMFCs  for  the  base  case,  k. 
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Fig.  7-Fig.  9.  Temperature  for  the  cross-section  of  Z= 0.05  m,  Z= 0.025  m,  Z= 0  m  for  the  single  phase  case,  k. 


The  lowest  temperature  could  be  observed  at  the  flow  inlet  where 
the  flows  with  constant  temperature  (353  K)  come  in,  and  the 
temperature  increases  slightly  along  the  positive  flow  direction  in 
the  cathode  channel  and  anode  gas  channel,  but  the  temperature 
changes  in  anode  channel  is  more  obvious  than  that  in  cath¬ 
ode  channel,  this  is  due  to  the  difference  of  thermal  conductivity 
between  air  and  hydrogen  (the  values  of  thermal  conductivity  for  air 
and  hydrogen  are  0.0242  and  0.16  respectively).  The  heat  generation 
is  due  to  the  electrochemical  reaction  and  Ohmic  heating.  There¬ 
fore,  the  highest  temperature  concentrates  in  the  catalyst  layers  in 
which  the  reactions  take  place  and  the  current  is  generated.  Espe¬ 
cially  in  the  area  neighboring  to  the  gas  channel,  the  temperature  is 
relative  high  for  the  weakness  of  the  gas  thermal  conductivity.  Due 
to  the  constant  temperature  for  the  coolant  boundary,  the  change  of 
temperature  distribution  through  the  current  collector  from  inlet 
to  outlet  is  small  enough  to  be  neglected  and  the  total  tempera¬ 
ture  difference  from  outlet  to  inlet  is  also  very  small.  To  investigate 
the  effect  of  water  saturation  on  the  temperature  distribution,  a 
single  phase  simulation  is  also  performed  in  the  PEMFCs  Module 
in  Fuent®  software.  Figs.  7-9  show  the  temperature  on  the  cross- 
section  of  Z= 0.05  m  (cathode  inlet),  Z= 0.025  m  (middle),  Z=0m 
(cathode  outlet)  for  the  single  phase  simulation.  It  can  be  seen  that 


the  relative  high  temperature  area  is  in  the  middle  cross-section 
which  is  reasonable  for  the  counter  flow  pattern  in  PEMFCs.  For 
demonstrating  the  effect  of  water  saturation,  the  temperature  in  the 
same  position  of  the  present  model  is  shown  in  Fig.  10  (Z= 0.05  m), 
Fig.  11  (Z= 0.025  m)  and  Fig.  12  (Z=  0). 

First,  the  comparison  between  the  cathode  inlet  for  the  two  cases 
(present  model  and  single  phase  model)  are  executed.  It  is  obvi¬ 
ous  that  the  temperature  distribution  are  nearly  the  same,  which 
is  consistent  with  the  fact  that  there  is  nearly  no  liquid  water  due 
to  the  vaporization  of  the  produced  water  and  little  accumulation. 
And  it  also  demonstrates  the  correctness  the  present  model  Sec¬ 
ondly,  the  temperature  difference  between  middle  cross-section  is 
checked,  and  it  is  found  that  the  relative  high  temperature  area  for 
the  present  two  phase  flow  model  is  greatly  larger  than  that  of  sin¬ 
gle  phase  simulation.  That  means  the  liquid  water  affects  the  heat 
transfer  in  the  PEMFCs,  and  it  hinders  the  heat  transfer  due  to  its  rel¬ 
ative  high  thermal  capacity  and  low  thermal  conductivity.  The  same 
difference  is  also  obtained  for  the  cross-section  ofZ=0,  where  the 
difference  is  more  obvious  for  the  increase  of  water  saturation.  It 
can  be  seen  that  the  relative  high  temperature  area  is  in  the  GDL  and 
catalyst  layer  neighboring  to  the  gas  channel,  but  not  neighboring 
to  the  current  collector  ribs,  although  where  the  water  saturation 
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Fig.  10-Fig.  12.  Temperature  for  the  cross-section  of  Z=  0.05  m,  Z= 0.025  m,  Z  =  0  m  for  the  base  case,  k. 
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is  high.  This  may  be  caused  by  the  fact  that  the  thermal  transfer 
capacity  of  cathode  gas  is  largely  small  than  that  of  solid  materials 
of  the  current  collector,  and  also  the  current  collector  boundary  is 
assumed  to  be  ideally  cooled  So,  although  the  water  saturation  is 
low  in  the  area  of  GDL  near  gas  channel  compared  with  that  of  near 
current  collector,  it  has  important  effect  on  the  heat  transfer,  which 
has  been  shown  in  Fig.  12. 

3.4.  The  effect  of  flow  pattern  on  the  temperature 

According  to  the  results  of  the  present  work  (Fig.  10-12),  It  can 
be  observed  that  the  maximum  hot  area  appears  at  the  middle  of 
the  PEMFCs  for  the  inlet  of  the  cathode  and  anode  are  both  kept  at 
the  constant  temperature  (353  K)  in  this  study  for  counter  flow  pat¬ 
tern.  To  make  insight  into  the  heat  transfer  difference  between  the 
counter  flow  pattern  and  coflow  pattern,  the  simulation  for  coflow 
pattern  is  also  performed  with  the  present  model.  Figs.  13-15  shows 
the  temperature  for  the  cross-section  of  Z=  0.05  m,  Z=  0.025  m  and 
Z=0.  And  Fig.  16  shows  the  slices  views  of  water  saturation  in  the 
PEMFCs.  It  can  be  seen  that  the  inlet  temperature  for  coflow  pat¬ 


tern  is  relative  low  than  that  of  counter  flow  pattern  (cathode  inlet), 
which  is  due  to  the  cooling  effect  of  the  inlet  gases  of  both  cath¬ 
ode  side  and  anode  side.  As  for  the  temperature  distribution  in 
the  middle  cross-section,  it  is  obvious  that  the  relative  high  tem¬ 
perature  area  of  counter  flow  pattern  is  large  than  that  of  coflow 
pattern.  And  for  the  outlet  cross-section,  the  maximum  tempera¬ 
ture  is  about  356.12  K  for  coflow  pattern,  which  is  larger  than  that 
of  counter  flow  pattern,  and  consistent  with  the  water  saturation 
distribution  in  Fig.  16.  So,  the  characters  of  the  temperature  dis¬ 
tribution  for  the  coflow  pattern  is  different  from  the  counter  flow 
pattern,  and  the  larger  temperature  difference  as  well  as  water  sat¬ 
uration  are  produced,  which  is  very  important  and  disadvantageous 
for  the  operation  of  PEMFCs. 

3.5.  The  effect  of  droplet  size  on  the  water  saturation 

In  the  two  phase  flow  for  the  liquid  droplets  or  solid  particles 
in  the  gas  phase,  the  size  of  the  particles  (droplets)  is  a  key  factor, 
which  directly  determine  the  two  phase  flow  pattern  and  disperse 
phase  (droplet  or  particle)  behavior.  It  can  be  seen  that  the  droplet 
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Fig.  13-Fig.  15.  Temperature  for  the  cross-section  of  Z= 0.05  m,  Z= 0.025  m,  Z= 0  m  for  the  coflow  pattern  case,  k. 
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Fig.  16.  Water  saturation  for  coflow  pattern. 


relax  time  increases  with  the  increase  of  the  droplet  size  according 
to  Eq.  ( 11 ),  which  means  the  large  droplet  needs  more  time  to  obtain 
the  same  velocity  with  the  small  droplet,  and,  if  the  droplet  is  large 
enough,  then  the  value  of  relax  time  is  very  big,  that  means  the 
dominated  motion  of  the  droplet  will  not  be  flowing  with  the  gas, 
so,  large  droplet  (the  ratio  between  the  characterize  size  and  the  gas 
channel  height  above  0.5)  [5]  for  the  very  low  gas  velocity  cannot  be 
considered  in  this  model,  the  same  for  the  other  two-fluid  model. 
And,  in  this  case,  VOF  model  is  more  suitable  for  the  two  phase 
simulation.  Figs.  17  and  18  show  (100  pm  and  200  pan,  the  droplet 
size  for  the  base  case  is  about  145  fxm)  the  water  saturation  for 
the  different  droplet  size.  It  is  obvious  that  the  water  saturation 
increases  with  the  increase  of  droplet  size,  which  means  the  large 
droplet  is  difficult  to  be  removed  and  has  more  effect  on  the  water 
saturation  in  GDL  and  catalyst  layer. 

However,  in  the  PEMFCs,  the  effect  of  the  gas  channel  proper¬ 
ties  on  the  liquid  removal  is  also  very  important,  which  has  been 
demonstrated  by  the  experiment,  so,  to  truly  simulate  the  two 
phase  behavior  in  PEMFCs,  the  parameters  reflecting  the  gas  chan¬ 
nel  properties  should  be  considered,  such  as  the  contact  angle  of  the 
gas  channel  wall,  but  all  of  the  models  for  the  whole  PEMFCs  in  the 
reference  fail  to  do  that.  In  the  future  work,  we  plan  to  explore  the 


198 


G.  He  et  al.  /  Journal  of  Power  Sources  194  (2009)  190-198 


used  multiphase  mixture  model  in  the  reference,  and  the  effect  of 
gas  diffusion  layer  (GDL)  properties  and  the  gas  drag  function  on 
the  liquid  water  removal  can  be  integrated  into  the  droplet  size 
as  a  parameter.  And  it  is  calculation  cost  effective  for  the  present 
model  compared  with  the  two-fluid  model  (which  also  includes  the 
effect  of  the  droplet  size).  The  three-dimensional  two  phase  sim¬ 
ulation  of  the  PEMFCs  with  a  straight  flow  field  is  performed,  and 
the  heat  transfer  is  also  included.  The  effect  of  droplet  size  on  the 
liquid  remove,  the  effect  of  liquid  water  on  the  heat  transfer  and  the 
effect  of  gas  flow  pattern  on  the  temperature  distribution  are  mainly 
investigated.  The  simulation  results  show  that  the  large  droplet  is 
hard  to  be  dragged  by  the  gas,  so  it  produces  large  water  satura¬ 
tion.  The  results  of  the  heat  transfer  show  that  the  liquid  water 
hinders  the  heat  transfer  in  the  GDL  and  catalyst  layer,  so  it  pro¬ 
duces  the  large  relative  high  temperature  area,  and  the  liquid  water 
has  important  effect  on  the  heat  transfer  in  the  area  of  GDL  and  cat¬ 
alyst  layer  near  the  channel  due  to  the  low  capacity  of  heat  transfer 
of  the  gas.  The  results  show  there  is  large  temperature  difference 
and  water  saturation  in  the  PEMFCs  operated  with  coflow  pattern 
compared  with  counter  flow  pattern,  which  is  disadvantageous  for 
the  operation  of  PEMFCs. 


Fig.  17.  Water  saturation  for  the  droplet  size  of  100  (Jim. 
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Fig.  18.  Water  saturation  for  the  droplet  size  of  200  p,m. 


relationship  between  the  gas  channel  contact  angle  and  the  size  of 
water  droplet  (or  film),  and  by  this  way,  the  effect  of  the  gas  channel 
can  be  included  in  the  present  model,  which  is  more  closely  to  the 
true  condition  of  the  PEMFCs. 


4.  Conclusions 

A  droplet  size  dependent  multiphase  mixture  model  is  devel¬ 
oped  for  the  two  phase  simulation  of  the  PEMFCs  in  this  paper,  and 
the  droplet  size  in  the  gas  channel  can  be  considered  as  a  param¬ 
eter  in  this  mixture  model  which  is  different  from  the  commonly 
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